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ABSTRACT 

We investigate the density-shear instability in Hall-MHD via numerical simulation of the full 
non-linear problem, in the context of magnetar activity. We confirm the development of the 
instability of a plane-parallel magnetic field with an appropriate intensity and electron den¬ 
sity profile, in accordance with analytic theory. We find that the instability also appears for a 
monotonically decreasing electron number density and magnetic field, a plane-parallel ana¬ 
logue of an azimuthal or meridional magnetic field in the crust of a magnetar. The growth rate 
of the instability depends on the Hall properties of the field (magnetic field intensity, electron 
number density and the corresponding scale-heights), while being insensitive to weak resistiv¬ 
ity. Since the Hall effect is the driving process for the evolution of the crustal magnetic field of 
magnetars, we argue that this instability is critical for systems containing strong meridional or 
azimuthal fields. We find that this process mediates the formation of localised structures with 
much stronger magnetic field than the average, which can lead to magnetar activity and accel¬ 
erate the dissipation of the field and consequently the production of Ohmic heating. Assuming 
a 5 x 10 14 G magnetic field at the base of crust, we anticipate that magnetic field as strong as 
10 15 G will easily develop in regions of typical size of a few 10 2 meters, containing magnetic 
energy of 10 43 erg, sufficient to power magnetar bursts. These active regions are more likely 
to appear in the magnetic equator where the tangential magnetic field is stronger. 
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1 INTRODUCTION 


The magnetic field evolution in the crust of neutron stars (NSs), in 
the magnetar regime, is mediated primarily by the Hall effect and 
Ohmic dissipation ( Jon es||1988[ |Goldreich & Reisenegger||1992| ). 
While the familiar picture of the Hall effect is that of the creation 
of a voltage across an electrical conductor, when a magnetic field 
is administered perpendicular to the current ( |Hall||1880| , NS ap¬ 
plications require that the feedback of the electric current onto the 
magnetic field is accounted for, leading to the realms of the Hall- 
Magnetohydrodynamics (Hall-MHD) description. 

As Hall-MHD provides a kinematic description, which does 
not correspond to an energy minimisation principle ( [Lyutikov] 
|2013| ), there is an ongoing debate regarding the stability of mag¬ 
netic configurations, turbulent cascade and the overall evolution¬ 
ary behaviour of a magnetic field in this context. Motivated by 
the mathematical similarity of the Hall-MHD equations with the 
vorticity equation in fluid dynamics, which is known to initiate 
turbulent cascade, it has been argued (Goldreich & Reisenegger 


|1992||B~iskamp et al.|1996||Cho & Lazarian|2009| > that a magnetic 

field should undergo Hall-induced turbulence. Plane parallel and 3- 
D cartesian box simulations ( [Warning & Hollerbach|[2009[ |2010) 
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demonstrated that while the magnetic field adopts a characteris¬ 
tic power spectrum once it evolves under Hall-MHD, its tempo¬ 
ral evolution in real space consists of frozen-in structures, whose 
time average is non-zero unlike normal turbulence. This result is in 
line with the consensus of axially-symmetric spherical-shell simu¬ 
lations of the magnetic field evolution in NS crusts where the Hall 
effect operates, with subdominant Ohmic dissipation. These simu- 


lations ([Hollerbach & Ru diger 2002 2004 ^ Pons et al.|2009]|Ko 


jima & Kisaka||2012| |Vigano et al.||2013| |Gourgouliatos & Cum 

' - b March *- ' ' . 


min g|2014a|bj|Marchant et al.|2014) " find that while the magnetic 

field may change drastically compared to its initial state as a re¬ 
sult of the Hall effect, the evolution saturates in a short time and the 
system relaxes to a particular spatial structure, a result recently con¬ 
firmed through 3-D spherical shell simulations ( [Wood & Holler- 1 
|bach|2015] >. 


In parallel to the turbulent cascade discussion, the question of 
Hall instability has been addressed. Although the Hall effect con¬ 
serves magnetic energy, it can potentially drive instability by trans¬ 
ferring energy from a Hall equilibrium state to a weaker perturbing 
field. Numerical and analytical studies have explored Hall insta¬ 
bility ([Rheinhardt & Geppert||2002[ |Rheinhardt et al.||2004[ |Pons| 
|& Geppert|2010| > in cartesian geometry. Quite remarkably, axially 

symmetric simulations in spherical shells did not find any evidence 
for the operation of Hall instability, i.e. when a state of Hall equilib- 
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rium ^Gourgouliatos et al.|2013} is chosen as an initial condition on 
an axially symmetric Hall simulation the system evolves because of 
Ohmic decay rather than the Hall effect ( [Marchant et al.|20l4) . 

| Wood et ah] ( |2014| ) studied analytically the density-shear in¬ 
stability for a unidirectional magnetic field. In this instability it is 
critical that both the magnetic field and the electron number density 
have strong gradients in the direction normal to the magnetic field. 
NS crusts are excellent environments for this instability to operate. 
They have a thickness of ~ 1km and host magnetic fields that could 
reach strengths of ~ 10 15 G for magnetars. The density at the base 
of the crust approaches the nuclear density ~ 10 14 g cnT 3 , and the 
Hall effect operates down to ~ 10 10 g cm -3 , as below this value, 
the effect of Lorentz forces becomes comparable to the breaking 
strain of the crust invalidating the Hall approach { |Gourgouliatos &| 
|Cumm ing 2015). As the chemical composition of the crust changes 
with depth and consequently the electron number fraction, the elec¬ 
tron number density in which the magnetic field evolved because of 
Hall-MHD ranges between ~ 10 36 cm -3 and 10 34 cm -3 (Cumming 
|et al.[2004| ). 

Transient activity of magnetars, in the form of bursts has 
been attributed to Hall evolution (Thompson & Murray|2001), via 
crust yielding. Elaborating on this scenario, Perna & Pons| ( |20lT] ) 
used axially symmetric Hall simulations to compare the magnetic 
stresses exerted on the crust to the breaking strain. They found that 
for initial poloidal fields B p - 8 x 10 14 G and toroidal B t = 2x 10 15 G 
magnetar activity is feasible, however a weaker initial magnetic 
field combination (B p = 2 x 10 14 G and B t = 10 15 G) leads only to 
sporadic bursts. Given that a substantial fraction of magn etars , have 
poloidal magnetic fields well below 5 x 10 14 G (Olausen & Kaspi 


|2014| ), it puts in question the validity of this scenario, given that 
magnetar behaviour has been observed by NSs with modest spin- 
down inferred dipole magnetic fields ( |Gavriil & Kaspi|2002| |Rea| 
|et al.| 2010, Sch olz et al.|2012| ). A possible solution to this puzzle is 
the presence of localised stronger magnetic fields compared to the 
large scale ones, a scenario that has been supported observationally 
fTiengo et al.|2013] >. In this work we show that the density-shear 
instability can severely deform the large scale structure of the mag¬ 
netic field in the crust of a NS and increase its intensity in areas of 
characteristic length-scale of a few 10 2 m, concentrating 10 43 erg of 
magnetic energy in the corresponding volume. This stronger mag¬ 
netic field exerts stresses in the crust that can lead to yielding and 
eventually to magnetar bursts. 

The plan of the paper is as follows: In Section 2 we provide the 
mathematical formulation of the problem. In Section 3 we present 
the numerical scheme and the initial conditions chosen. We discuss 
the results in Section 4. We consider applications to magnetar ac¬ 
tivity in Section 5. We conclude in Section 6. 


2 MATHEMATICAL FORMULATION 

In the electron-MHD limit of the Hall effect, the electron fluid ve¬ 
locity v e is related to the electric current density by j = 
where n e is the electron number density, c and e are the speed of 
light and the electron elementary charge. Then, from Ampere’s law, 
the electric current density is j - x B, where B is the mag¬ 
netic induction. We can safely neglect Maxwell’s correction as the 
velocities involved are non-relativistic. The electron velocity be¬ 
comes v e = V x B. Assuming some finite conductivity cr, the 
electric field reads E = -v e x B/c + j/cr. Finally we substitute into 


Faraday’s law to obtain the induction equation: 


dB 

dt 


-Vx 


—!— (V x B) x B + V x B 
4nen e 4 ncr 


( 1 ) 


The first term in the right hand side of equation {!]) describes the 
evolution of the magnetic field under the influence of the Hall ef¬ 
fect, while the second term describes Ohmic dissipation. We define 
a timescale for the Hall effect t H = 4yr ^ L , where L is the typical 

length scale of the problem, while for Ohmic decay it is t 0 = \ 

the ratio of t 0 lt H gives the dimensionless Magnetic Reynolds num¬ 
ber R b - —, also referred to as the Hall Parameter. 

Having assumed a plane-parallel geometry, the system is in¬ 
variant to translations in the y direction and the quantities depend 
only on x and z. We then express the magnetic field in terms of two 
scalar functions: 


B = VW(x, z)xy + B y (x, z)y , 


( 2 ) 


which is by construction divergence free. Substituting expression 
{2} into the induction equation {T}, we obtain two coupled differen¬ 
tial equations for B y and ¥: 


d^ 

dt 

dB y 

dt 


■ (V5 v x v) • VY + ——V 2v P, 

A y J ) dnrr 


4nn e Q ' y ' 4ncr 

4;re \ \ n e 


X y j • + By (Vn~ ] x y) • YB V 


X-^B y -a-^B y .Vo-). 


(3) 


(4) 


We switch to dimensionless quantities, keeping the same notation. 


dt [\ 


+ R 


-l 

B 


(VByXy)- V'F + flg'V 2 '?, 
(' V 2 B y - o-~ l WB y ■ Vo-) . 


xy)-VBy 


(5) 

( 6 ) 


Appropriate profiles of n e (x ) are imposed; cr is taken to be a con¬ 
stant. In our simulation the unit time is tn, while the Ohmic dis¬ 
sipation time is R B tH- Because of the varying electron density and 
magnetic field throughout the domain, Hall evolution may develop 
substantially faster than this timescale. 


3 NUMERICAL SIMULATION 

We integrate the full non-linear equations ^ and |6), using Euler’s 
method, in a uniform grid i, k so that x = idx and z = kdz. We 
apply periodic boundary conditions in z. Regarding the x boundary 
condition we use two setups. First, we assume a vacuum in either 
side of the x boundary, by fitting a current-free magnetic field for 
x < -1 and x > 1 (BC1); in this setup the large scale magnetic 
field and the electron density profile are symmetric about the axis 
x — 0. This is used to confirm the occurrence of the instability and 
to compare with the analytical model. In the second setup, we use 
the vacuum boundary condition for x > 1, while for x < -1 we as¬ 
sume that there is no magnetic field penetrating that boundary, by 
setting ^ - By - 0atx = -l (BC2). This condition is more restric¬ 
tive than the Meissner superconductor boundary condition regard¬ 
ing the By component ( [Hollerbach & Rudiger |2004| , nevertheless, 
it is a good approximation once R B » 1, as is the case here. In this 
configuration, the electron number density has its maximum values 
at x = -1 and decrease monotonically, resembling the structure of a 
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NS crust. We implement these boundary conditions by using an ap¬ 
propriate set of ghost points. We have tested the results in different 
resolution levels to ensure their validity. We use a courant condi¬ 
tion that adjusts the timestep depending on the maximum electron 
velocity. 

It has been shown analytically that the density-shear instabil¬ 
ity occurs when the (dimensionless) magnetic field and the elec¬ 
tron number density profiles are chosen so that B z (x) = n e (x) = 
sech y (x), where y is some positive constant ( Wood et al. 2014). 
Taking the asymptotic limit for y —> 0 for the expression sech r (x/y) 
and y —» oo of the expression sech r (x/y 1/2 ) we find respectively the 
backgrounds B - n 0 = exp(-|jc|) and B - n - exp(-x 2 /2). In this 
work we focus in the gaussian profile because of its smoothness. 
We have also run some simulations using the absolute value profile 
to validate the occurrence of the instability. 

We implement these profiles as follows. The absolute value 
profile where the initial condition for the magnetic field is 
B = B 0 (exp(-|.x/L 5 |) + e B ) Z and n Q = n 0 (exp(-|x/LJ) + e n ), 
and the gaussian profile with initial magnetic field B - 
B 0 (2n~ 1/2 exp(-x 2 /L 2 ) + e 5 ) z and n e = n 0 (exp {-x 2 /L 2 ) + e n ). We 
superimpose a perturbation term b = -6b(cos(k p z)x + sin (k p z)y). 
We have included a uniform background field e B B 0 z and a uniform 
background density e n n 0 , with e„ « 1 and 6b e B B (h to ensure 
that the perturbing magnetic field and currents are always subdom¬ 
inant compared to the background field; the typical values used for 
the background field is 10 ~ 2 B 0 and the perturbation 10 ~ 4 B 0 . The 
above profiles are used with the boundary condition BC1 where the 
system is symmetric with respect to x = 0. We also used a trans¬ 
lated version of the gaussian profile where x —» x + 1, imposing 
BC2, in this case the code dissipates some energy to force the per¬ 
turbation on B y to satisfy the boundary condition at x = - 1 . 

We also run simulations using a pseudo-spectral parallel code. 
This code implements the second order, Runge-Kutta ETD time¬ 
stepping scheme described in Cox & Matthews (2002), and has 
been modified to integrate the Hall-MHD equations. The main dif¬ 
ference is that we employ periodic boundary conditions both in x 
and z boundaries, unlike the grid-based one which assumes vacuum 
or the non-penetrating field condition in the x direction. 


4 RESULTS 

We have explored various combinations of the parameters. A sum¬ 
mary is shown in Table 1, where we provide information on the 
initial conditions, and the resulting instability. In addition to the 
quantities already defined we give the wavenumber of the fastest 
growing mode k t , the corresponding growth timescale r and the 
resolution used. 

We have confirmed that a uniform magnetic field on a gaussian 
density background (and vice versa, runs SI and S2) does not lead 
to any unstable mode. We have run simulations using the gaussian 
profile and BC1 boundary conditions for a broad combination of 
parameters (G), we considered no background uniform field e B = 0, 
while keeping the other quantities the same (G9) which also gave 
rise to the instability. We have also used a smaller number of simu¬ 
lations using the absolute value profile (A) and BC1 boundary con¬ 
ditions. Applying BC2, we run two simulations (C) using the trans¬ 
lated gaussian profile B - B 0 (2n 1/2 exp(-(v + 1) 2 /L 2 ) + e B ^j z and 

n Q = n 0 (exp(-(* + 1) 2 /L 2 ) + e n ). 

We confirm the development of the instability once the scale- 
height of the magnetic field and the density variation are compa¬ 
rable within a range of a few, and the resistivity is weak. The sys¬ 


tem undergoes some adjustment, followed by exponential growth 
of the instability, see Fig. |T| Once the instability fully develops, its 
energy content is comparable to that of the background magnetic 
field, with the overall structure being deformed, Fig. [2] At this point 
large electron velocities develop, leading to a very small timestep 
forcing us to stop our calculation, a numerical limitation known to 
exist in explicit Eulerian Hall-MHD simulations ( |Falle|2003| . The 
wavenumber of the fastest growing mode of the instability depends 
on the scale height of the magnetic field and electron number den¬ 
sity, being inversely proportional to them once L B - L n , see for 
instance the kf s of GO and G7. However if the scale-heights of the 
magnetic field and the density are not equal, the evolution becomes 
more complex, with the magnetic field needing extra time to adjust 
to the density background before the instability starts growing (GO 
vs G6). The growth rate is proportional to the strength of the mag¬ 
netic field, i.e. GO, G1 and G2, where the ratio of the respective r’s 
is 0.5 while Bq is increasing by 2. In the limit of strong resistivity 
(i.e. G4), the instability may be suppressed, without dominating the 
overall evolution, even though there is some modest growth at the 
beginning. 


Repeating the analysis of |Wood et al.|j20T4) , for the gaus¬ 
sian profile under our normalisation for L B - L n - L, we find that 
the growth rate co 2 = B 2 k 2 (2 - L 2 k 2 )/(L 2 /ra 0 ), with the maximum 
rate occurring for k t = L~ x giving oj max = B 0 /(L 2 y/ra 0 ). As the 
wavenumber of the fastest growing mode of the instability is small, 
it is affected by the size of the simulation box, being forced to be a 
multiple of n because of the periodic boundary conditions imposed. 
Having assumed L n - L B = 0.1, n 0 = 1 = B 0 (G8), the simulation 
gives ki = 3n = 9.42 versus an analytical value of k t - 0.1 -1 = 10 
and a corresponding growth timescale r(= cd~ 1 ) = 0.0138 versus 
an analytical prediction of 0.0178. This deviation is due to nu¬ 
merical constraints and also to the superimposed uniform magnetic 
field and background density which are not present in the analyti¬ 
cal model. Because of numerical limitations we have not been able 
to set a strong constraint on the maximum and minimum ratio of 
L B /L n where the instability appears, except for the fact that there 
is no instability for uniform magnetic field or density. To investi¬ 
gate that, it would require either L B «: L n or L B » L n and both 
of them to be much smaller than the size of the box, leading to a 
calculation that ranges over a few orders of magnitude. In the sim¬ 
ulations where we used the absolute value profile, even at a very 
low resistivity (Al) the system undergoes some significant decay 
as the currents are very strong around x = 0, which slows down the 
growth of the instability, while choices of higher resistivity (A2) 
prevent its development entirely. 


Similar behaviour is found when BC2 is applied, Fig.[3] Given 
that there is a rigid boundary at x = -1 there is a significant growth 
of the B x and B z components of the magnetic field, because of the 
compression of the magnetic field lines against the boundary, com¬ 
pared to the other case, where the main effect of the instability was 
to kink the structure of the field. These results are in broad qual¬ 
itative agreement with the linear calculation of |Rheinhardt et alT] 
( 2004), done in a similar setup. Using the parallel code, we simu¬ 
lated the development of the instability for R B = 50 on a collocation 
grid with 256 2 points and the same initial conditions and density 
profiles as the GO run (see run SP in Table 1), finding the same 
behaviour. 
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Figure 1. The ratio of the energy in the B y component over the total mag¬ 
netic energy for some characteristic runs (for details on the parameters refer 
to Table 1). With the exception of the highly dissipative model G4, all other 
models undergo some rearrangement of the perturbing field which is is fol¬ 
lowed by exponential growth of the instability. 


Table 1. Simulations summary. The S runs have either a uniform mag¬ 
netic field or density background, the G runs utilise the gaussian profile, the 
A runs the absolute value profile while the C runs utilise a monotinically 
decreasing magnetic field and electron number density profile with BC2 
boundary conditions. 


NAME 

Bo 

Lb 

L n 

r b 

kp/n 

kiln 

t / 10 -2 

Resol. 

SI 

1 

0.1 

- 

0 

10 

- 

- 

200 2 

S2 

1 

- 

0.1 

0 

10 

- 

- 

200 2 

GO 

1 

0.1 

0.1 

0.01 

5 

3 

1.16 

200 2 

G1 

2 

0.1 

0.1 

0.01 

5 

3 

0.613 

200 2 

G2 

4 

0.1 

0.1 

0.01 

5 

3 

0.323 

200 2 

G3 

1 

0.1 

0.1 

0.05 

5 

2 

1.36 

200 2 

G4 

1 

0.1 

0.1 

0.1 

5 

1 

- 

200 2 

G5 

1 

0.1 

0.05 

0.001 

20 

4 

0.183 

200 2 

G6 

1 

0.05 

0.1 

0.001 

10 

2 

2.76 

200 2 

G7 

1 

0.05 

0.05 

0.001 

20 

6 

0.282 

200 2 

G8 

1 

0.1 

0.1 

0 

5 

3 

1.38 

200 2 

G9 

1 

0.1 

0.1 

0 

10 

3 

1.36 

200 2 

Al 

1 

0.1 

0.1 

0.001 

10 

2 

3.44 

100 2 

A2 

1 

0.1 

0.1 

0.005 

10 

1 

- 

100 2 

Cl 

1 

0.1 

0.1 

0 

10 

2 

1.76 

200 2 

C2 

1 

0.15 

0.1 

0 

10 

2 

0.731 

100 2 

SP 

1 

0.1 

0.1 

0.02 

10 

3 

1.83 

256 2 




- 1.0 - 0.5 0.0 0.5 1.0 

x 

Figure 2. The structure of the magnetic field at the beginning of the simu¬ 
lation (top panel) and once the instability has fully developed at t - 0.15^// 
(bottom panel), for the simulation GO. The B x and B z components are plot¬ 
ted in black, while the B y component is shown in colour. 


Figure 3. The structure of the magnetic field once the instability has fully 
developed at t - 0.175, for the simulation Cl, where BC2 is used. 

5 APPLICATION TO MAGNETAR ACTIVITY 

Consider a magnetar, whose magnetic field and electron number 
density at the base of the crust are 5 x 10 14 G, 10 36 cm -3 respectively 
and crust thickness is 1km. Using this normalisation in eqns ([5} 
and & we find that the unit time of the simulation corresponds to 
~ 10 5 years in the NS’s life. The scale-height for the electron num¬ 
ber density is ~ 0.1km, thus we expect a growth time scale of ~ 10 3 
years, as we have found that r ~ 10 -2 . Even under the conservative 
assumption of a magnetic field with a tiny amount of energy being 
in the perturbing small scale field, it is only a mater of a few thou¬ 
sand years for this instability to create a strong localised magnetic 
field, exceeding the intensity of the background field by a factor 
of 2 and giving rise to magnetic fields ~ 10 15 G for this setup. The 
size of these structures is comparable to the scale-height multiplied 
by 2/r, thus they are expected to be ~ 0.6km, each one of them 
containing magnetic energy 10 43 erg, which is sufficient to power 
magnetar busts. This effect is caused entirely by the redistribution 
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of the magnetic field via the Hall effect without appealing to the 
generation of any extra magnetic flux. Even if a moderately strong 
large scale magnetic field is present (5 x 10 13 - 10 14 G), this instabil¬ 
ity leads to the formation of pockets of magnetic field significantly 
exceeding the average value. In our simulations we found that these 
features typically develop near the base of the crust rather than the 
surface, as in our initial condition we have chosen an exponentially 
decreasing profile. While it is possible that such features may de¬ 
velop closer to the surface, this is a question to be answered con¬ 
clusively by future more realistic simulations. This is particularly 
interesting in the context of recent observations of strong localised 
magnetic features such as the one observed in SGR 0418+5729 
ffiengo et al.|2013f and the 0.2-0.7 km hotspot implied by surface 
emission modelling in the same system ( [Guillot et al.|2015| >. 

As the components of the magnetic field which are parallel to 
layers of constant density are susceptible to this instability we ex¬ 
pect the non-radial magnetic field (meridional and toroidal) to con¬ 
tribute the most. In a typical large-scale poloidal dipole magnetic 
field structure, the meridional component is stronger away from 
the poles, making these instabilities more likely to develop in mid¬ 
latitudes and in the equatorial region, with respect to the magnetic 
dipole axis. Thus we expect bursts triggered through this mecha¬ 
nism to provide energy away from the poles, leading to the appear¬ 
ance of hot spots in the form of subpulses, of the same frequency 
yet different phase compared to the main pulse which is likely to 
be associated to the magnetic pole. Recent observations show that 
bursts are evenly distributed in spin phase ( [Collazzi et al.|2015). In 
any case, the complexity of heat transport within the crust ( [Brown 
& Cumming 2009) and the size of the active region (Baubock et al. 
2015)7 are critical for the observational appearance of these fea¬ 
tures. 


6 CONCLUSIONS 

In this work we have confirmed numerically the development of the 
density-shear instability in a plane-parallel geometry. In particular, 
we have found that the instability appears when the scale height 
of the magnetic field and the electron number density are compa¬ 
rable, with the growth timescale depending on the intensity of the 
magnetic field, the electron number density and the relevant scale- 
heights. This instability also appears in a monotonically decreasing 
electron number density and magnetic field, a structure that encap¬ 
sulates the basic characteristics of a NS crust. We conclude that 
the density-shear instability can lead to the formation of localised 
strong magnetic fields, with the typical size of these areas being a 
few times the scale height. Realistic NS studies need to go beyond 
this plane-parallel geometry, test the appearance of this instability 
in a 3-D calculation and investigate its evolution after it has fully 
developed, a task which is out the capacity of the current numerical 
scheme. Nevertheless, it is likely that a natural NS configuration 
can host an appropriate magnetic field geometry that will give rise 
to this instability and provide an efficient mechanism for powering 
magnetar activity with a weaker overall magnetic field. 
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